function [C_est, t_est, cost, flag] = dls_pnp(p, z)
% DLS-PnP:
%
% This function performs the DLS-PnP method introduced at ICCV 2011
% Joel A. Hesch and Stergios I. Roumeliotis. "A direct least-squares (dls) 
% solution for PnP". In Proc. of the Int. Conf. on Computer Vision, 
% Barcelona, Spain, November 6-13, 2011. 
%
% inputs: 
% p: 3xN vector of 3D known point features
% z: 2xN vector of correpsonding image measurements (calibrated)

% Check the inputs
if size(z,1) > size(z,2) || size(p,1) > size(p,2)
    fprintf('Usage: dls_pnp(p,z) \n p: 3xN matrix of 3D points \n z: 2xN matrix of corresponding 2D image measurements (normalized pixel coordinates)')
end

% make z into unit vectors from normalized pixel coords
z = [z; ones(1,size(z,2))];
z = z./ repmat(sqrt(sum(z.*z,1)),3,1);

% some preliminaries
flag = 0;
N = size(z,2);
 
% build coeff matrix
% An intermediate matrix, the inverse of what is called "H" in the paper
% (see eq. 25)
H = zeros(3);
for i = 1:N
    H = H + eye(3) - z(:,i)*z(:,i)';
end

A = zeros(3,9);
for i = 1:N
   A = A + (z(:,i)*z(:,i)' - eye(3)) * LeftMultVec(p(:,i));
end
A = H\A;

D = zeros(9);
for i = 1:N
    D = D + (LeftMultVec(p(:,i)) + A)' * (eye(3) - z(:,i)*z(:,i)') * (LeftMultVec(p(:,i)) + A);
end

f1coeff = [2*D(1,6) - 2*D(1,8) + 2*D(5,6) - 2*D(5,8) + 2*D(6,1) + 2*D(6,5) + 2*D(6,9) - 2*D(8,1) - 2*D(8,5) - 2*D(8,9) + 2*D(9,6) - 2*D(9,8); % constant term
          (6*D(1,2) + 6*D(1,4) + 6*D(2,1) - 6*D(2,5) - 6*D(2,9) + 6*D(4,1) - 6*D(4,5) - 6*D(4,9) - 6*D(5,2) - 6*D(5,4) - 6*D(9,2) - 6*D(9,4)); % s1^2  * s2
          (4*D(1,7) - 4*D(1,3) + 8*D(2,6) - 8*D(2,8) - 4*D(3,1) + 4*D(3,5) + 4*D(3,9) + 8*D(4,6) - 8*D(4,8) + 4*D(5,3) - 4*D(5,7) + 8*D(6,2) + 8*D(6,4) + 4*D(7,1) - 4*D(7,5) - 4*D(7,9) - 8*D(8,2) - 8*D(8,4) + 4*D(9,3) - 4*D(9,7)); % s1 * s2
          (4*D(1,2) - 4*D(1,4) + 4*D(2,1) - 4*D(2,5) - 4*D(2,9) + 8*D(3,6) - 8*D(3,8) - 4*D(4,1) + 4*D(4,5) + 4*D(4,9) - 4*D(5,2) + 4*D(5,4) + 8*D(6,3) + 8*D(6,7) + 8*D(7,6) - 8*D(7,8) - 8*D(8,3) - 8*D(8,7) - 4*D(9,2) + 4*D(9,4)); % s1 * s3
          (8*D(2,2) - 8*D(3,3) - 8*D(4,4) + 8*D(6,6) + 8*D(7,7) - 8*D(8,8)); % s2 * s3
          (4*D(2,6) - 2*D(1,7) - 2*D(1,3) + 4*D(2,8) - 2*D(3,1) + 2*D(3,5) - 2*D(3,9) + 4*D(4,6) + 4*D(4,8) + 2*D(5,3) + 2*D(5,7) + 4*D(6,2) + 4*D(6,4) - 2*D(7,1) + 2*D(7,5) - 2*D(7,9) + 4*D(8,2) + 4*D(8,4) - 2*D(9,3) - 2*D(9,7)); % s2^2 * s3
          (2*D(2,5) - 2*D(1,4) - 2*D(2,1) - 2*D(1,2) - 2*D(2,9) - 2*D(4,1) + 2*D(4,5) - 2*D(4,9) + 2*D(5,2) + 2*D(5,4) - 2*D(9,2) - 2*D(9,4)); % s2^3
          (4*D(1,9) - 4*D(1,1) + 8*D(3,3) + 8*D(3,7) + 4*D(5,5) + 8*D(7,3) + 8*D(7,7) + 4*D(9,1) - 4*D(9,9)); % s1 * s3^2
          (4*D(1,1) - 4*D(5,5) - 4*D(5,9) + 8*D(6,6) - 8*D(6,8) - 8*D(8,6) + 8*D(8,8) - 4*D(9,5) - 4*D(9,9)); % s1
          (2*D(1,3) + 2*D(1,7) + 4*D(2,6) - 4*D(2,8) + 2*D(3,1) + 2*D(3,5) + 2*D(3,9) - 4*D(4,6) + 4*D(4,8) + 2*D(5,3) + 2*D(5,7) + 4*D(6,2) - 4*D(6,4) + 2*D(7,1) + 2*D(7,5) + 2*D(7,9) - 4*D(8,2) + 4*D(8,4) + 2*D(9,3) + 2*D(9,7)); % s3
          (2*D(1,2) + 2*D(1,4) + 2*D(2,1) + 2*D(2,5) + 2*D(2,9) - 4*D(3,6) + 4*D(3,8) + 2*D(4,1) + 2*D(4,5) + 2*D(4,9) + 2*D(5,2) + 2*D(5,4) - 4*D(6,3) + 4*D(6,7) + 4*D(7,6) - 4*D(7,8) + 4*D(8,3) - 4*D(8,7) + 2*D(9,2) + 2*D(9,4)); % s2
          (2*D(2,9) - 2*D(1,4) - 2*D(2,1) - 2*D(2,5) - 2*D(1,2) + 4*D(3,6) + 4*D(3,8) - 2*D(4,1) - 2*D(4,5) + 2*D(4,9) - 2*D(5,2) - 2*D(5,4) + 4*D(6,3) + 4*D(6,7) + 4*D(7,6) + 4*D(7,8) + 4*D(8,3) + 4*D(8,7) + 2*D(9,2) + 2*D(9,4)); % s2 * s3^2
          (6*D(1,6) - 6*D(1,8) - 6*D(5,6) + 6*D(5,8) + 6*D(6,1) - 6*D(6,5) - 6*D(6,9) - 6*D(8,1) + 6*D(8,5) + 6*D(8,9) - 6*D(9,6) + 6*D(9,8)); % s1^2
          (2*D(1,8) - 2*D(1,6) + 4*D(2,3) + 4*D(2,7) + 4*D(3,2) - 4*D(3,4) - 4*D(4,3) - 4*D(4,7) - 2*D(5,6) + 2*D(5,8) - 2*D(6,1) - 2*D(6,5) + 2*D(6,9) + 4*D(7,2) - 4*D(7,4) + 2*D(8,1) + 2*D(8,5) - 2*D(8,9) + 2*D(9,6) - 2*D(9,8)); % s3^2
          (2*D(1,8) - 2*D(1,6) - 4*D(2,3) + 4*D(2,7) - 4*D(3,2) - 4*D(3,4) - 4*D(4,3) + 4*D(4,7) + 2*D(5,6) - 2*D(5,8) - 2*D(6,1) + 2*D(6,5) - 2*D(6,9) + 4*D(7,2) + 4*D(7,4) + 2*D(8,1) - 2*D(8,5) + 2*D(8,9) - 2*D(9,6) + 2*D(9,8)); % s2^2
          (2*D(3,9) - 2*D(1,7) - 2*D(3,1) - 2*D(3,5) - 2*D(1,3) - 2*D(5,3) - 2*D(5,7) - 2*D(7,1) - 2*D(7,5) + 2*D(7,9) + 2*D(9,3) + 2*D(9,7)); % s3^3
          (4*D(1,6) + 4*D(1,8) + 8*D(2,3) + 8*D(2,7) + 8*D(3,2) + 8*D(3,4) + 8*D(4,3) + 8*D(4,7) - 4*D(5,6) - 4*D(5,8) + 4*D(6,1) - 4*D(6,5) - 4*D(6,9) + 8*D(7,2) + 8*D(7,4) + 4*D(8,1) - 4*D(8,5) - 4*D(8,9) - 4*D(9,6) - 4*D(9,8)); % s1 * s2 * s3
          (4*D(1,5) - 4*D(1,1) + 8*D(2,2) + 8*D(2,4) + 8*D(4,2) + 8*D(4,4) + 4*D(5,1) - 4*D(5,5) + 4*D(9,9)); % s1 * s2^2
          (6*D(1,3) + 6*D(1,7) + 6*D(3,1) - 6*D(3,5) - 6*D(3,9) - 6*D(5,3) - 6*D(5,7) + 6*D(7,1) - 6*D(7,5) - 6*D(7,9) - 6*D(9,3) - 6*D(9,7)); % s1^2 * s3
          (4*D(1,1) - 4*D(1,5) - 4*D(1,9) - 4*D(5,1) + 4*D(5,5) + 4*D(5,9) - 4*D(9,1) + 4*D(9,5) + 4*D(9,9))]; % s1^3

f2coeff = [- 2*D(1,3) + 2*D(1,7) - 2*D(3,1) - 2*D(3,5) - 2*D(3,9) - 2*D(5,3) + 2*D(5,7) + 2*D(7,1) + 2*D(7,5) + 2*D(7,9) - 2*D(9,3) + 2*D(9,7); % constant term
            (4*D(1,5) - 4*D(1,1) + 8*D(2,2) + 8*D(2,4) + 8*D(4,2) + 8*D(4,4) + 4*D(5,1) - 4*D(5,5) + 4*D(9,9)); % s1^2  * s2
            (4*D(1,8) - 4*D(1,6) - 8*D(2,3) + 8*D(2,7) - 8*D(3,2) - 8*D(3,4) - 8*D(4,3) + 8*D(4,7) + 4*D(5,6) - 4*D(5,8) - 4*D(6,1) + 4*D(6,5) - 4*D(6,9) + 8*D(7,2) + 8*D(7,4) + 4*D(8,1) - 4*D(8,5) + 4*D(8,9) - 4*D(9,6) + 4*D(9,8)); % s1 * s2
            (8*D(2,2) - 8*D(3,3) - 8*D(4,4) + 8*D(6,6) + 8*D(7,7) - 8*D(8,8)); % s1 * s3
            (4*D(1,4) - 4*D(1,2) - 4*D(2,1) + 4*D(2,5) - 4*D(2,9) - 8*D(3,6) - 8*D(3,8) + 4*D(4,1) - 4*D(4,5) + 4*D(4,9) + 4*D(5,2) - 4*D(5,4) - 8*D(6,3) + 8*D(6,7) + 8*D(7,6) + 8*D(7,8) - 8*D(8,3) + 8*D(8,7) - 4*D(9,2) + 4*D(9,4)); % s2 * s3
            (6*D(5,6) - 6*D(1,8) - 6*D(1,6) + 6*D(5,8) - 6*D(6,1) + 6*D(6,5) - 6*D(6,9) - 6*D(8,1) + 6*D(8,5) - 6*D(8,9) - 6*D(9,6) - 6*D(9,8)); % s2^2 * s3
            (4*D(1,1) - 4*D(1,5) + 4*D(1,9) - 4*D(5,1) + 4*D(5,5) - 4*D(5,9) + 4*D(9,1) - 4*D(9,5) + 4*D(9,9)); % s2^3
            (2*D(2,9) - 2*D(1,4) - 2*D(2,1) - 2*D(2,5) - 2*D(1,2) + 4*D(3,6) + 4*D(3,8) - 2*D(4,1) - 2*D(4,5) + 2*D(4,9) - 2*D(5,2) - 2*D(5,4) + 4*D(6,3) + 4*D(6,7) + 4*D(7,6) + 4*D(7,8) + 4*D(8,3) + 4*D(8,7) + 2*D(9,2) + 2*D(9,4)); % s1 * s3^2
            (2*D(1,2) + 2*D(1,4) + 2*D(2,1) + 2*D(2,5) + 2*D(2,9) - 4*D(3,6) + 4*D(3,8) + 2*D(4,1) + 2*D(4,5) + 2*D(4,9) + 2*D(5,2) + 2*D(5,4) - 4*D(6,3) + 4*D(6,7) + 4*D(7,6) - 4*D(7,8) + 4*D(8,3) - 4*D(8,7) + 2*D(9,2) + 2*D(9,4)); % s1
            (2*D(1,6) + 2*D(1,8) - 4*D(2,3) + 4*D(2,7) - 4*D(3,2) + 4*D(3,4) + 4*D(4,3) - 4*D(4,7) + 2*D(5,6) + 2*D(5,8) + 2*D(6,1) + 2*D(6,5) + 2*D(6,9) + 4*D(7,2) - 4*D(7,4) + 2*D(8,1) + 2*D(8,5) + 2*D(8,9) + 2*D(9,6) + 2*D(9,8)); % s3
            (8*D(3,3) - 4*D(1,9) - 4*D(1,1) - 8*D(3,7) + 4*D(5,5) - 8*D(7,3) + 8*D(7,7) - 4*D(9,1) - 4*D(9,9)); % s2
            (4*D(1,1) - 4*D(5,5) + 4*D(5,9) + 8*D(6,6) + 8*D(6,8) + 8*D(8,6) + 8*D(8,8) + 4*D(9,5) - 4*D(9,9)); % s2 * s3^2
            (2*D(1,7) - 2*D(1,3) + 4*D(2,6) - 4*D(2,8) - 2*D(3,1) + 2*D(3,5) + 2*D(3,9) + 4*D(4,6) - 4*D(4,8) + 2*D(5,3) - 2*D(5,7) + 4*D(6,2) + 4*D(6,4) + 2*D(7,1) - 2*D(7,5) - 2*D(7,9) - 4*D(8,2) - 4*D(8,4) + 2*D(9,3) - 2*D(9,7)); % s1^2
            (2*D(1,3) - 2*D(1,7) + 4*D(2,6) + 4*D(2,8) + 2*D(3,1) + 2*D(3,5) - 2*D(3,9) - 4*D(4,6) - 4*D(4,8) + 2*D(5,3) - 2*D(5,7) + 4*D(6,2) - 4*D(6,4) - 2*D(7,1) - 2*D(7,5) + 2*D(7,9) + 4*D(8,2) - 4*D(8,4) - 2*D(9,3) + 2*D(9,7)); % s3^2
            (6*D(1,3) - 6*D(1,7) + 6*D(3,1) - 6*D(3,5) + 6*D(3,9) - 6*D(5,3) + 6*D(5,7) - 6*D(7,1) + 6*D(7,5) - 6*D(7,9) + 6*D(9,3) - 6*D(9,7)); % s2^2
            (2*D(6,9) - 2*D(1,8) - 2*D(5,6) - 2*D(5,8) - 2*D(6,1) - 2*D(6,5) - 2*D(1,6) - 2*D(8,1) - 2*D(8,5) + 2*D(8,9) + 2*D(9,6) + 2*D(9,8)); % s3^3
            (8*D(2,6) - 4*D(1,7) - 4*D(1,3) + 8*D(2,8) - 4*D(3,1) + 4*D(3,5) - 4*D(3,9) + 8*D(4,6) + 8*D(4,8) + 4*D(5,3) + 4*D(5,7) + 8*D(6,2) + 8*D(6,4) - 4*D(7,1) + 4*D(7,5) - 4*D(7,9) + 8*D(8,2) + 8*D(8,4) - 4*D(9,3) - 4*D(9,7)); % s1 * s2 * s3
            (6*D(2,5) - 6*D(1,4) - 6*D(2,1) - 6*D(1,2) - 6*D(2,9) - 6*D(4,1) + 6*D(4,5) - 6*D(4,9) + 6*D(5,2) + 6*D(5,4) - 6*D(9,2) - 6*D(9,4)); % s1 * s2^2
            (2*D(1,6) + 2*D(1,8) + 4*D(2,3) + 4*D(2,7) + 4*D(3,2) + 4*D(3,4) + 4*D(4,3) + 4*D(4,7) - 2*D(5,6) - 2*D(5,8) + 2*D(6,1) - 2*D(6,5) - 2*D(6,9) + 4*D(7,2) + 4*D(7,4) + 2*D(8,1) - 2*D(8,5) - 2*D(8,9) - 2*D(9,6) - 2*D(9,8)); % s1^2 * s3
            (2*D(1,2) + 2*D(1,4) + 2*D(2,1) - 2*D(2,5) - 2*D(2,9) + 2*D(4,1) - 2*D(4,5) - 2*D(4,9) - 2*D(5,2) - 2*D(5,4) - 2*D(9,2) - 2*D(9,4))]; % s1^3

f3coeff = [2*D(1,2) - 2*D(1,4) + 2*D(2,1) + 2*D(2,5) + 2*D(2,9) - 2*D(4,1) - 2*D(4,5) - 2*D(4,9) + 2*D(5,2) - 2*D(5,4) + 2*D(9,2) - 2*D(9,4); % constant term
          (2*D(1,6) + 2*D(1,8) + 4*D(2,3) + 4*D(2,7) + 4*D(3,2) + 4*D(3,4) + 4*D(4,3) + 4*D(4,7) - 2*D(5,6) - 2*D(5,8) + 2*D(6,1) - 2*D(6,5) - 2*D(6,9) + 4*D(7,2) + 4*D(7,4) + 2*D(8,1) - 2*D(8,5) - 2*D(8,9) - 2*D(9,6) - 2*D(9,8)); % s1^2  * s2
          (8*D(2,2) - 8*D(3,3) - 8*D(4,4) + 8*D(6,6) + 8*D(7,7) - 8*D(8,8)); % s1 * s2
          (4*D(1,8) - 4*D(1,6) + 8*D(2,3) + 8*D(2,7) + 8*D(3,2) - 8*D(3,4) - 8*D(4,3) - 8*D(4,7) - 4*D(5,6) + 4*D(5,8) - 4*D(6,1) - 4*D(6,5) + 4*D(6,9) + 8*D(7,2) - 8*D(7,4) + 4*D(8,1) + 4*D(8,5) - 4*D(8,9) + 4*D(9,6) - 4*D(9,8)); % s1 * s3
          (4*D(1,3) - 4*D(1,7) + 8*D(2,6) + 8*D(2,8) + 4*D(3,1) + 4*D(3,5) - 4*D(3,9) - 8*D(4,6) - 8*D(4,8) + 4*D(5,3) - 4*D(5,7) + 8*D(6,2) - 8*D(6,4) - 4*D(7,1) - 4*D(7,5) + 4*D(7,9) + 8*D(8,2) - 8*D(8,4) - 4*D(9,3) + 4*D(9,7)); % s2 * s3
          (4*D(1,1) - 4*D(5,5) + 4*D(5,9) + 8*D(6,6) + 8*D(6,8) + 8*D(8,6) + 8*D(8,8) + 4*D(9,5) - 4*D(9,9)); % s2^2 * s3
          (2*D(5,6) - 2*D(1,8) - 2*D(1,6) + 2*D(5,8) - 2*D(6,1) + 2*D(6,5) - 2*D(6,9) - 2*D(8,1) + 2*D(8,5) - 2*D(8,9) - 2*D(9,6) - 2*D(9,8)); % s2^3
          (6*D(3,9) - 6*D(1,7) - 6*D(3,1) - 6*D(3,5) - 6*D(1,3) - 6*D(5,3) - 6*D(5,7) - 6*D(7,1) - 6*D(7,5) + 6*D(7,9) + 6*D(9,3) + 6*D(9,7)); % s1 * s3^2
          (2*D(1,3) + 2*D(1,7) + 4*D(2,6) - 4*D(2,8) + 2*D(3,1) + 2*D(3,5) + 2*D(3,9) - 4*D(4,6) + 4*D(4,8) + 2*D(5,3) + 2*D(5,7) + 4*D(6,2) - 4*D(6,4) + 2*D(7,1) + 2*D(7,5) + 2*D(7,9) - 4*D(8,2) + 4*D(8,4) + 2*D(9,3) + 2*D(9,7)); % s1
          (8*D(2,2) - 4*D(1,5) - 4*D(1,1) - 8*D(2,4) - 8*D(4,2) + 8*D(4,4) - 4*D(5,1) - 4*D(5,5) + 4*D(9,9)); % s3
          (2*D(1,6) + 2*D(1,8) - 4*D(2,3) + 4*D(2,7) - 4*D(3,2) + 4*D(3,4) + 4*D(4,3) - 4*D(4,7) + 2*D(5,6) + 2*D(5,8) + 2*D(6,1) + 2*D(6,5) + 2*D(6,9) + 4*D(7,2) - 4*D(7,4) + 2*D(8,1) + 2*D(8,5) + 2*D(8,9) + 2*D(9,6) + 2*D(9,8)); % s2
          (6*D(6,9) - 6*D(1,8) - 6*D(5,6) - 6*D(5,8) - 6*D(6,1) - 6*D(6,5) - 6*D(1,6) - 6*D(8,1) - 6*D(8,5) + 6*D(8,9) + 6*D(9,6) + 6*D(9,8)); % s2 * s3^2
          (2*D(1,2) - 2*D(1,4) + 2*D(2,1) - 2*D(2,5) - 2*D(2,9) + 4*D(3,6) - 4*D(3,8) - 2*D(4,1) + 2*D(4,5) + 2*D(4,9) - 2*D(5,2) + 2*D(5,4) + 4*D(6,3) + 4*D(6,7) + 4*D(7,6) - 4*D(7,8) - 4*D(8,3) - 4*D(8,7) - 2*D(9,2) + 2*D(9,4)); % s1^2
          (6*D(1,4) - 6*D(1,2) - 6*D(2,1) - 6*D(2,5) + 6*D(2,9) + 6*D(4,1) + 6*D(4,5) - 6*D(4,9) - 6*D(5,2) + 6*D(5,4) + 6*D(9,2) - 6*D(9,4)); % s3^2
          (2*D(1,4) - 2*D(1,2) - 2*D(2,1) + 2*D(2,5) - 2*D(2,9) - 4*D(3,6) - 4*D(3,8) + 2*D(4,1) - 2*D(4,5) + 2*D(4,9) + 2*D(5,2) - 2*D(5,4) - 4*D(6,3) + 4*D(6,7) + 4*D(7,6) + 4*D(7,8) - 4*D(8,3) + 4*D(8,7) - 2*D(9,2) + 2*D(9,4)); % s2^2
          (4*D(1,1) + 4*D(1,5) - 4*D(1,9) + 4*D(5,1) + 4*D(5,5) - 4*D(5,9) - 4*D(9,1) - 4*D(9,5) + 4*D(9,9)); % s3^3
          (4*D(2,9) - 4*D(1,4) - 4*D(2,1) - 4*D(2,5) - 4*D(1,2) + 8*D(3,6) + 8*D(3,8) - 4*D(4,1) - 4*D(4,5) + 4*D(4,9) - 4*D(5,2) - 4*D(5,4) + 8*D(6,3) + 8*D(6,7) + 8*D(7,6) + 8*D(7,8) + 8*D(8,3) + 8*D(8,7) + 4*D(9,2) + 4*D(9,4)); % s1 * s2 * s3
          (4*D(2,6) - 2*D(1,7) - 2*D(1,3) + 4*D(2,8) - 2*D(3,1) + 2*D(3,5) - 2*D(3,9) + 4*D(4,6) + 4*D(4,8) + 2*D(5,3) + 2*D(5,7) + 4*D(6,2) + 4*D(6,4) - 2*D(7,1) + 2*D(7,5) - 2*D(7,9) + 4*D(8,2) + 4*D(8,4) - 2*D(9,3) - 2*D(9,7)); % s1 * s2^2
          (4*D(1,9) - 4*D(1,1) + 8*D(3,3) + 8*D(3,7) + 4*D(5,5) + 8*D(7,3) + 8*D(7,7) + 4*D(9,1) - 4*D(9,9)); % s1^2 * s3
          (2*D(1,3) + 2*D(1,7) + 2*D(3,1) - 2*D(3,5) - 2*D(3,9) - 2*D(5,3) - 2*D(5,7) + 2*D(7,1) - 2*D(7,5) - 2*D(7,9) - 2*D(9,3) - 2*D(9,7))]; % s1^3
      
% Construct the Macaulay matrix
% u0 = round(randn(1)*100);
% u1 = round(randn(1)*100);
% u2 = round(randn(1)*100);
% u3 = round(randn(1)*100);
%M2 = cayley_LS_M(f1coeff, f2coeff, f3coeff,u0,u1,u2,u3);

u = round(randn(4,1) * 100);
M2 = cayley_LS_M(f1coeff, f2coeff, f3coeff,u);

% construct the multiplication matrix via schur compliment of the Macaulay
% matrix
Mtilde = M2(1:27,1:27) - M2(1:27,28:120)/M2(28:120,28:120)*M2(28:120,1:27);

[V,~] = eig(Mtilde); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Now check the solutions                                               %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% extract the optimal solutions from the eigen decomposition of the
% Multiplication matrix
sols = zeros(3,27);
cost = zeros(27,1);
i = 1;
for k = 1:27
    V(:,k) = V(:,k)/V(1,k);
    if (imag(V(2,k)) == 0)
        
        stmp = V([10 4 2],k);

        H = Hessian(f1coeff, f2coeff, f3coeff, stmp);
        if eig(H) > 0
            sols(:,i) = stmp;
            
            Cbar = cayley2rotbar(stmp);  
            CbarVec = Cbar';
            CbarVec = CbarVec(:);
            cost(i) = CbarVec' * D * CbarVec;
            i = i+1;
        end
    end
end
sols = sols(:,1:i-1);
cost = cost(1:i-1);

C_est = zeros(3,3,size(sols,2));
t_est = zeros(3, size(sols,2));

for j = 1:size(sols,2)

    % recover the optimal orientation
    C_est(:,:,j) = 1/(1 + sols(:,j)' * sols(:,j)) * cayley2rotbar(sols(:,j));
    
    A2 = zeros(3);
    for i = 1:N
        A2 = A2 + eye(3) - z(:,i)*z(:,i)';
    end
    
    b2 = zeros(3,1);
    
    for i = 1:N
        b2 = b2 + (z(:,i)*z(:,i)' - eye(3)) * C_est(:,:,j) * p(:,i);
    end
    
    % recover the optimal translation
    t_est(:,j) = A2\b2;
end

% check that the points are infront of the center of perspectivity
sols_valid = [];
for k = 1:size(sols,2)
    cam_points = C_est(:,:,k) * p + repmat(t_est(:,k),1, length(p));
    
    if isempty(find(cam_points(3,:) < 0))
        sols_valid = [sols_valid; k];
    end
    
end
t_est = t_est(:,sols_valid);
C_est = C_est(:,:,sols_valid);
cost = cost(sols_valid);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Some helper functions                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function C = cayley2rotbar(s)

C = ( (1-s'*s) * eye(3) + 2 * skewsymm(s) + 2 * (s * s'))';

end

function C = skewsymm(X1)
% generates skew symmetric matrix

C = [0      , -X1(3) ,  X1(2)
    X1(3) , 0      , -X1(1)
    -X1(2) , X1(1)  ,     0];
end

function M = LeftMultVec(v)
% R * p = LeftMultVec(p) * vec(R)
M = [v' zeros(1,6);
     zeros(1,3) v' zeros(1,3);
     zeros(1,6) v'];
end

function M = cayley_LS_M(a,b,c,u) %,u1,u2,u3)
% Construct the Macaulay resultant matrix

M = [u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) 0; u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) a(10) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(10) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(10) 0; 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) a(14) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 b(10) 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(14) 0 0 0 0 0 c(1) 0 0 0 0 0 0 0 0 0 c(10) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(14) 0; u(3) 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 0 0 b(11) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) c(1); 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) a(5) 0 0 0 0 0 0 0 a(1) 0 0 0 0 0 a(10) 0 0 0 0 b(11) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 b(10) 0 0 0 0 0 0 b(5) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) c(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(5) c(10); 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(5) a(12) 0 0 0 0 0 a(1) 0 a(10) 0 0 0 0 0 a(14) 0 a(11) 0 0 b(5) 0 0 0 0 0 0 0 b(1) 0 0 b(11) 0 0 0 0 0 b(10) 0 0 0 0 b(14) 0 0 0 0 0 0 b(12) 0 0 0 0 0 c(11) 0 0 0 0 0 0 0 0 0 c(5) c(10) 0 0 0 0 0 0 0 0 0 0 c(1) 0 0 0 0 0 0 c(12) c(14); 0 0 0 u(3) 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 0 0 a(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) b(11) 0 0 0 0 0 0 b(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) 0 0 0 0 0 0 0 0 0 0 c(15) c(11); 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) 0 0 0 a(15) a(6) 0 0 0 0 0 0 0 a(11) 0 a(1) 0 0 0 a(5) 0 0 0 0 b(15) 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 0 0 0 b(11) 0 0 0 b(10) b(5) 0 0 0 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(15) c(11) 0 0 0 0 0 0 c(10) 0 0 0 0 c(1) 0 0 0 0 0 c(6) c(5); 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) 0 0 0 a(6) 0 0 0 0 0 0 a(11) 0 a(5) 0 a(10) a(1) 0 0 a(12) 0 a(15) 0 0 b(6) 0 0 0 0 0 0 0 b(11) b(10) 0 b(15) b(1) 0 0 0 0 b(5) 0 0 0 b(14) b(12) 0 0 0 0 0 0 0 0 0 0 0 0 c(15) 0 0 0 0 0 0 0 0 0 c(6) c(5) 0 c(1) 0 0 0 0 c(14) 0 0 0 c(11) c(10) 0 0 0 0 0 0 c(12); u(2) 0 0 0 0 0 0 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(9) a(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) b(1) 0 0 0 c(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(9) 0; 0 u(2) 0 0 0 0 0 0 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) a(9) a(4) a(10) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) 0 0 0 0 b(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(4) b(10) 0 0 0 c(10) 0 0 0 0 0 0 0 0 0 0 c(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) c(4) 0; 0 0 u(2) 0 0 0 0 0 0 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 0 0 a(10) a(4) a(8) a(14) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(9) 0 0 b(4) 0 0 b(1) 0 b(10) 0 0 0 0 0 b(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(8) b(14) 0 0 0 c(14) c(9) 0 0 0 0 0 0 0 0 0 c(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) 0 0 c(10) c(8) 0; 0 0 0 u(2) 0 0 0 0 0 u(3) 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(3) a(11) 0 0 a(1) 0 0 0 0 0 0 0 0 0 a(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) 0 0 b(1) 0 0 0 b(3) b(11) 0 0 0 c(11) 0 0 0 0 0 0 0 c(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(3) c(9); 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) a(3) a(17) a(5) 0 0 a(10) 0 0 0 a(9) 0 0 0 a(1) 0 a(4) 0 0 0 0 b(3) 0 0 0 0 b(11) b(1) 0 0 0 0 0 0 0 0 0 0 b(9) 0 0 0 0 b(4) 0 0 b(10) 0 0 0 b(17) b(5) 0 0 0 c(5) 0 c(1) 0 0 0 0 0 c(10) 0 0 c(3) c(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) c(17) c(4); 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) 0 0 0 0 a(5) a(17) 0 a(12) 0 0 a(14) 0 a(9) a(1) a(4) 0 0 0 a(10) 0 a(8) 0 a(3) 0 0 b(17) 0 b(1) b(11) 0 b(5) b(10) 0 b(9) 0 0 b(3) 0 0 0 0 0 b(4) 0 0 0 0 b(8) 0 0 b(14) 0 0 0 0 b(12) 0 0 0 c(12) c(3) c(10) 0 0 0 0 0 c(14) 0 0 c(17) c(4) 0 0 0 0 0 c(1) 0 0 0 0 c(9) 0 0 c(11) 0 0 c(5) 0 c(8); 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) a(9) 0 0 0 0 a(18) a(15) 0 0 a(11) 0 0 0 0 0 0 0 0 0 a(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) b(9) b(3) 0 0 b(11) 0 0 0 b(18) b(15) 0 0 0 c(15) 0 0 0 0 0 c(1) 0 c(11) 0 0 0 0 0 0 0 0 0 0 c(9) 0 0 0 0 0 0 0 0 0 0 c(18) c(3); 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 a(10) a(4) 0 0 a(15) a(18) 0 a(6) 0 0 a(5) 0 0 0 a(3) a(1) a(9) 0 a(11) 0 a(17) 0 0 0 0 b(18) 0 0 0 0 b(15) b(11) 0 0 b(9) 0 0 0 0 b(1) 0 0 b(3) 0 0 b(10) b(4) b(17) 0 0 b(5) 0 0 0 0 b(6) 0 0 0 c(6) 0 c(11) 0 0 0 c(10) 0 c(5) c(1) 0 c(18) c(3) 0 0 0 0 0 0 c(4) 0 0 0 0 c(9) 0 0 0 0 c(15) 0 c(17); 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 a(15) a(14) a(8) 0 0 a(6) 0 0 0 0 0 a(12) 0 a(3) a(11) a(17) a(10) a(4) a(9) a(5) 0 0 0 a(18) 0 0 0 0 b(11) b(15) 0 b(6) b(5) 0 b(3) b(4) 0 b(18) b(9) 0 b(10) 0 0 b(17) 0 0 b(14) b(8) 0 0 0 b(12) 0 0 0 0 0 0 0 0 0 c(18) c(5) 0 0 0 c(14) 0 c(12) c(10) 0 0 c(17) 0 c(9) 0 0 0 c(11) c(8) 0 0 0 c(3) c(4) 0 c(15) 0 0 c(6) 0 0; 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(13) a(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) b(13) b(9) 0 0 c(1) c(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(13) 0; 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 0 0 0 a(1) a(9) a(13) a(19) a(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) b(13) 0 0 0 0 b(9) 0 b(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(10) b(19) b(4) 0 0 c(10) c(4) 0 0 0 0 0 0 0 0 0 0 c(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) c(9) c(19) 0; 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 u(4) u(1) 0 0 0 0 0 0 a(1) a(9) 0 0 0 a(10) a(4) a(19) 0 a(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(13) 0 a(14) b(19) b(1) 0 b(9) 0 b(4) 0 b(10) 0 0 0 b(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(14) 0 b(8) 0 0 c(14) c(8) c(13) 0 0 0 0 0 0 0 0 0 c(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) c(9) 0 c(10) c(4) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 0 u(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(2) a(3) 0 a(1) a(9) 0 0 0 0 0 0 0 0 0 a(13) 0 0 0 a(11) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(13) 0 b(1) b(9) 0 0 b(11) b(2) b(3) 0 0 c(11) c(3) 0 0 0 c(1) 0 0 0 c(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(2) c(13); 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 0 0 0 0 0 0 a(11) a(3) a(2) 0 a(17) 0 a(10) a(4) a(1) 0 0 a(13) 0 0 0 a(9) 0 a(19) 0 0 0 a(5) b(2) 0 0 0 0 b(3) b(9) b(11) 0 0 0 0 0 0 0 0 0 b(13) b(1) 0 0 0 b(19) 0 b(10) b(4) 0 0 b(5) 0 b(17) 0 0 c(5) c(17) 0 c(9) 0 c(10) 0 0 c(1) c(4) 0 0 c(2) c(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) c(3) 0 c(19); 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 a(11) a(3) 0 0 0 a(5) a(17) 0 0 0 0 a(14) a(8) a(10) a(13) a(9) a(19) 0 0 0 a(4) 0 0 0 a(2) 0 a(12) 0 b(11) b(9) b(3) 0 b(17) b(4) b(5) b(13) 0 0 b(2) 0 0 0 0 0 b(19) b(10) 0 0 0 0 0 b(14) b(8) 0 0 b(12) 0 0 0 0 c(12) 0 c(2) c(4) 0 c(14) 0 0 c(10) c(8) 0 0 0 c(19) 0 0 0 0 0 c(9) 0 0 0 0 c(13) 0 c(11) c(3) 0 c(5) c(17) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 0 u(1) 0 0 0 0 a(9) a(13) 0 0 0 0 0 a(18) 0 a(11) a(3) 0 0 0 0 0 0 0 0 0 a(2) 0 0 a(1) a(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) b(9) b(13) b(2) 0 b(11) b(3) 0 0 b(15) 0 b(18) 0 0 c(15) c(18) 0 0 0 c(11) c(1) c(9) 0 c(3) 0 0 0 0 0 0 0 0 0 0 c(13) 0 0 0 0 0 0 0 0 0 0 0 c(2); 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) 0 0 0 a(4) a(19) 0 a(15) a(18) 0 0 0 0 a(5) a(17) a(11) 0 0 a(2) a(9) a(13) 0 a(3) 0 0 0 0 a(10) a(6) 0 0 0 0 0 b(18) b(3) b(15) 0 b(13) 0 0 0 0 b(9) 0 0 b(2) b(11) b(10) b(4) b(19) 0 0 b(5) b(17) 0 0 b(6) 0 0 0 0 c(6) 0 0 c(3) 0 c(5) c(10) c(4) c(11) c(17) c(9) 0 0 c(2) 0 0 0 0 0 0 c(19) 0 0 0 0 c(13) 0 0 0 c(15) c(18) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 u(3) 0 u(4) u(1) a(15) a(18) a(8) 0 0 a(6) 0 0 0 0 0 a(12) 0 a(5) a(2) a(3) 0 a(4) a(19) a(13) a(17) 0 0 0 0 a(14) 0 0 b(15) b(3) b(18) 0 0 b(17) b(6) b(2) b(19) 0 0 b(13) 0 b(4) 0 0 0 b(5) b(14) b(8) 0 0 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 c(17) 0 c(12) c(14) c(8) c(5) 0 c(4) 0 0 0 0 c(13) 0 0 0 c(3) 0 0 0 0 c(2) c(19) c(15) c(18) 0 c(6) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) a(3) 0 0 0 0 0 a(7) 0 0 a(15) 0 0 0 0 0 0 0 0 0 a(18) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) b(1) 0 0 0 b(11) b(3) b(18) 0 0 b(15) 0 0 0 0 b(7) 0 0 0 c(7) 0 0 c(1) 0 0 c(11) 0 c(15) 0 0 0 0 0 0 0 0 0 0 c(3) 0 0 c(9) 0 0 0 0 0 0 0 0 c(18); 0 0 0 0 0 0 u(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(11) 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 0 b(11) b(15) 0 0 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) 0 0 c(1) 0 0 0 0 0 0 0 c(7) c(15); 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(3) 0 0 0 0 a(3) a(2) 0 0 0 0 0 0 0 a(15) a(18) 0 0 0 0 0 0 0 0 0 0 0 0 a(11) a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(13) b(9) 0 0 b(11) b(3) b(2) 0 0 b(15) b(18) 0 0 b(7) 0 0 0 0 c(7) 0 0 0 c(9) c(15) c(11) c(3) 0 c(18) 0 0 0 0 0 0 0 0 0 0 c(2) 0 0 c(13) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(18) 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 a(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) 0 0 b(2) 0 0 0 b(18) 0 0 0 b(7) 0 0 0 c(7) 0 0 0 0 0 c(20) 0 c(2) 0 0 0 0 c(18) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(2) 0 0 0 0 0 0 0 a(15) a(18) 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 a(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) b(13) 0 0 b(3) b(2) 0 0 b(15) b(18) 0 b(7) 0 0 0 0 0 c(7) 0 0 0 0 c(13) c(18) c(3) c(2) 0 0 0 c(15) 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(9) 0 a(13) 0 0 a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(13) b(9) b(20) 0 0 c(9) c(13) c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(15) a(18) 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) 0 0 0 0 b(3) b(11) 0 0 0 b(15) b(18) 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 c(11) 0 0 c(15) 0 c(7) 0 0 0 0 0 0 c(9) 0 0 0 c(18) 0 0 c(3) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(18) 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(15) 0 0 0 0 0 0 0 0 0 0 0 b(13) 0 0 0 0 b(2) b(3) 0 0 b(15) b(18) 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 c(3) c(7) c(15) c(18) 0 0 0 0 0 0 0 0 c(13) 0 0 0 0 0 0 c(2) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(18) 0 0 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 0 0 b(2) 0 0 b(18) 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 c(2) 0 c(18) 0 0 0 0 c(7) 0 0 0 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 u(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 0 0 0 0 0 a(10) 0 a(14) 0 0 0 0 0 a(16) 0 a(5) 0 0 b(12) 0 0 0 0 0 0 0 b(10) 0 0 b(5) 0 0 0 0 0 b(14) 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 c(5) 0 0 0 0 0 0 0 0 0 c(12) c(14) c(1) 0 0 0 0 0 0 0 c(11) 0 c(10) 0 0 0 0 0 0 0 c(16); 0 0 u(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) 0 0 b(14) 0 0 0 0 0 0 0 0 0 0 b(10) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 c(10) 0 0 0 0 0 0 0 0 0 c(14) 0 0 0 0 0 0 0 0 0 c(1) 0 0 0 0 0 0 0 0 c(16) 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 0 0 0 b(11) 0 0 0 0 0 b(15) b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) 0 0 0 c(15) 0 0 c(11) 0 0 0 0 0 0 0 0 c(7); 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(8) 0 0 0 0 0 b(14) 0 b(16) 0 0 0 0 0 b(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(9) 0 0 c(10) c(4) 0 0 0 0 c(14) 0 0 c(16) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) a(8) 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(14) 0 b(8) 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(13) 0 0 c(4) c(19) 0 0 0 c(14) c(8) 0 c(16) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(11) 0 0 0 0 b(15) 0 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) 0 0 0 c(7) 0 0 c(15) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(8) 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(8) 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 c(19) 0 0 0 0 c(8) 0 c(16) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(3) 0 0 0 0 b(18) b(15) 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(15) 0 0 c(7) 0 0 0 0 0 0 0 0 c(3) 0 0 0 0 0 0 c(18) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(4) 0 0 c(14) c(8) 0 0 0 0 c(16) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(19) 0 0 c(8) 0 0 0 0 c(16) 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 b(2) 0 0 0 0 0 b(18) 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(18) 0 c(7) 0 0 0 0 0 0 0 0 0 c(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 a(18) 0 0 0 0 0 0 0 a(6) 0 0 0 0 0 0 0 0 0 0 0 b(19) 0 0 b(2) b(18) 0 b(17) 0 b(7) b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(17) 0 c(6) 0 c(7) 0 c(18) 0 0 0 0 0 c(19) c(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(6) 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 a(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(15) b(10) 0 0 b(11) 0 b(5) 0 b(7) 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) 0 0 c(10) c(11) 0 0 c(6) 0 0 c(5) 0 c(15) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 0 0 a(16) a(14) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(12) b(16) 0 0 b(14) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(5) c(14) 0 0 c(15) 0 0 0 c(6) 0 c(12) c(16) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(14) 0 0 0 c(5) 0 0 0 c(12) 0 c(16) 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(18) 0 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) 0 0 0 0 0 0 0 0 0 0 0 c(18) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(6) 0 0 0 a(12) a(5) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(6) b(12) 0 0 b(5) b(14) 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(15) c(5) 0 c(14) 0 0 0 0 c(7) c(16) c(6) c(12) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(15) 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(15) 0 0 0 0 0 0 c(7) 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(14) 0 0 0 c(16) 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(17) 0 0 b(18) b(7) 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(6) 0 0 0 0 0 c(7) 0 0 0 0 0 c(17) c(18) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 0 c(16) 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(6) 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(6) c(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(12) 0 b(7) b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) c(12) c(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 u(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 a(5) 0 a(12) 0 a(14) a(10) 0 0 0 0 a(6) 0 0 0 0 0 0 0 0 0 0 b(5) b(14) 0 b(6) b(10) 0 0 0 0 b(12) 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 c(6) 0 0 0 0 0 0 0 0 0 0 c(12) c(11) c(10) 0 0 0 0 c(16) 0 c(15) 0 c(5) c(14) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 u(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 0 0 0 0 0 0 0 0 0 a(15) 0 a(6) 0 a(5) a(11) 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 b(15) b(5) 0 b(7) b(11) b(10) 0 b(14) 0 b(6) 0 0 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) 0 0 0 0 0 0 0 0 0 0 c(6) 0 c(11) 0 c(10) 0 0 c(12) 0 0 c(14) c(15) c(5) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 0 0 c(17) c(16) 0 c(12) 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(12) b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(6) c(12) 0 c(16) c(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 b(6) b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) c(6) c(16) c(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(12) c(16) 0 0 c(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 c(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) 0 a(16) 0 0 0 0 0 0 0 a(12) 0 0 0 0 0 0 0 0 0 0 b(14) 0 0 b(12) 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(12) 0 0 0 0 0 0 0 0 0 0 c(16) c(10) 0 0 0 c(11) 0 0 0 c(5) 0 c(14) 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(2) 0 0 0 0 0 0 0 0 0 0 a(18) 0 0 0 a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 b(2) 0 0 0 b(18) 0 0 0 c(18) 0 0 0 0 0 0 0 c(20) 0 0 0 0 c(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 a(9) a(13) 0 0 a(10) a(4) a(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) a(20) 0 a(8) 0 b(9) 0 b(13) b(10) b(19) 0 b(4) 0 0 0 b(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(14) 0 b(8) 0 0 0 c(14) c(8) 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(9) c(13) c(10) c(4) c(19) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(7) b(5) 0 0 b(15) 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(5) c(15) 0 0 0 0 0 c(6) 0 c(7) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 a(6) a(15) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(7) b(6) b(14) 0 b(15) b(5) 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(15) c(14) c(5) 0 0 0 0 0 c(12) c(7) c(6) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(13) 0 a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) b(13) 0 0 0 c(13) c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(3) 0 0 0 a(17) 0 0 a(7) 0 0 0 0 0 a(6) 0 a(15) 0 0 0 a(3) a(2) 0 a(18) 0 0 0 0 a(5) 0 0 0 0 0 0 0 b(18) b(7) 0 b(2) 0 0 0 b(13) b(3) b(19) b(4) 0 b(15) b(5) b(17) 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 c(18) c(4) c(6) c(5) c(17) c(15) 0 c(3) 0 0 0 0 0 0 c(13) 0 0 0 0 0 c(19) 0 c(2) 0 0 0 c(7) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(2) a(1) a(9) a(13) 0 0 0 0 0 0 0 0 0 a(20) a(11) 0 0 a(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) b(1) b(9) b(13) b(11) 0 b(3) 0 b(2) 0 c(11) c(3) c(2) 0 0 0 c(9) 0 0 0 c(13) 0 c(1) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20); 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 b(14) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(14) 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 0 0 c(1) 0 0 0 c(10) 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(17) 0 0 0 a(12) 0 0 0 0 0 a(16) 0 0 a(8) 0 a(19) 0 0 0 0 0 0 0 0 0 0 0 0 b(17) b(19) 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 b(8) 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 0 c(16) 0 0 c(20) 0 0 0 0 c(19) 0 c(2) 0 0 0 0 c(17) 0 c(12) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 a(16) 0 a(8) 0 0 0 0 0 0 0 0 0 b(12) 0 0 0 0 0 0 0 0 0 b(8) 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 c(17) c(8) 0 0 c(18) c(12) 0 c(6) 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(5) 0 0 0 0 a(12) 0 0 0 0 0 a(16) 0 a(4) a(10) a(8) 0 0 0 a(14) 0 0 0 a(17) 0 0 0 0 b(10) b(5) 0 b(12) b(14) 0 b(4) 0 0 b(17) 0 0 0 0 0 b(8) 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 c(17) c(14) 0 0 0 0 0 c(16) 0 0 0 c(8) c(9) 0 0 0 0 c(10) 0 c(11) c(3) 0 c(4) 0 0 c(5) 0 0 c(12) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(3) a(2) 0 0 0 0 a(4) a(19) 0 a(13) 0 0 0 0 0 0 a(20) a(5) 0 a(17) 0 0 0 0 0 0 0 b(3) 0 b(20) b(2) 0 0 0 0 0 0 0 0 0 0 b(13) 0 0 0 0 b(4) b(19) 0 b(17) b(5) 0 0 0 c(5) c(17) 0 0 0 c(20) 0 c(19) 0 0 c(13) 0 0 c(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(3) c(2) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(15) a(18) 0 a(7) 0 0 0 0 0 0 0 0 0 0 0 0 b(7) 0 0 b(18) b(4) 0 0 b(3) b(15) b(17) b(5) 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) c(5) 0 0 c(6) 0 0 c(15) 0 0 0 0 0 c(4) c(3) 0 0 0 0 0 c(17) 0 c(18) 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) a(18) 0 0 0 a(6) 0 0 0 0 0 a(12) 0 0 a(17) 0 a(2) 0 a(19) 0 a(20) 0 0 0 0 0 a(8) 0 0 b(18) b(2) 0 b(6) 0 0 0 0 0 0 0 b(20) 0 b(19) 0 0 0 b(17) b(8) 0 0 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 c(17) 0 c(19) c(12) 0 0 0 c(20) 0 0 0 c(2) 0 0 0 0 0 0 c(18) 0 c(6) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(3) 0 0 0 0 0 0 0 0 0 0 0 0 a(5) a(17) 0 0 a(7) 0 0 0 0 0 a(6) 0 0 0 a(18) a(11) a(3) 0 a(15) 0 0 0 0 0 0 0 0 0 0 0 b(7) b(15) 0 0 b(3) 0 0 0 b(9) b(11) b(4) b(10) b(18) 0 0 b(5) b(17) 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 0 c(15) c(10) 0 0 c(5) 0 c(6) c(11) 0 0 c(18) 0 0 0 c(9) 0 0 c(17) 0 0 c(4) 0 c(3) 0 0 0 0 c(7) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 a(19) 0 a(15) a(18) 0 0 0 0 a(5) a(17) 0 a(3) 0 0 0 a(13) a(20) 0 a(2) 0 0 a(6) 0 a(4) 0 0 0 0 0 b(15) 0 b(2) b(18) 0 b(20) 0 0 0 0 b(13) 0 0 0 b(3) b(4) b(19) 0 0 b(5) b(17) 0 b(6) 0 0 0 0 0 c(6) 0 0 0 c(2) 0 c(17) c(4) c(19) c(3) 0 c(13) c(5) 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 c(15) c(18) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(18) 0 0 0 0 0 a(17) 0 0 a(2) 0 0 0 a(20) 0 0 0 a(6) 0 0 0 a(19) 0 0 0 0 0 b(18) 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 b(2) b(19) 0 0 0 b(17) 0 0 0 b(6) 0 0 0 c(6) 0 0 0 0 0 0 0 c(19) 0 c(2) 0 c(20) c(17) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(18) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 0 a(11) a(3) a(2) 0 0 0 a(10) a(4) a(19) a(9) 0 0 a(20) 0 0 0 a(13) 0 0 a(5) 0 0 a(17) 0 0 0 0 b(11) b(2) b(13) b(3) 0 0 0 0 0 0 0 0 0 b(20) b(9) 0 0 0 0 b(10) b(4) b(19) b(5) 0 b(17) 0 0 0 c(5) c(17) 0 0 c(13) 0 c(4) 0 0 c(9) c(19) 0 c(10) 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(11) c(3) c(2) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(2) 0 0 0 a(17) 0 0 0 0 0 a(8) 0 0 a(19) 0 a(20) 0 0 0 0 0 a(12) 0 0 0 0 0 0 b(2) b(20) 0 b(17) 0 0 0 0 0 0 0 0 0 0 0 0 0 b(19) 0 0 0 0 b(8) 0 0 0 b(12) 0 0 0 c(12) 0 0 0 0 0 0 0 0 0 c(19) 0 0 c(8) 0 0 0 0 0 0 0 c(20) 0 0 0 0 0 0 c(2) 0 c(17) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 0 0 0 0 a(6) 0 0 a(18) 0 0 0 a(2) 0 0 0 0 0 0 0 a(17) 0 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 b(20) b(2) 0 b(19) 0 b(18) b(17) 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 c(19) 0 c(17) 0 c(18) 0 c(2) c(6) 0 0 0 0 0 c(20) 0 0 0 0 0 0 0 0 0 0 c(7) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 a(8) 0 0 0 0 0 0 0 0 0 0 0 0 b(12) b(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 0 0 c(19) 0 0 0 c(2) c(8) 0 c(17) 0 0 0 0 c(12) 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 a(3) a(2) 0 0 a(5) a(17) 0 0 0 0 a(14) a(8) 0 a(4) a(20) a(13) 0 0 0 0 a(19) 0 0 a(12) 0 0 0 0 b(3) b(13) b(2) b(5) 0 b(19) b(17) b(20) 0 0 0 0 0 0 0 0 0 b(4) 0 0 0 0 b(14) b(8) 0 b(12) 0 0 0 0 0 c(12) 0 0 0 c(19) 0 c(8) 0 0 c(4) 0 0 c(14) 0 0 0 0 0 0 0 c(13) 0 0 0 0 c(20) 0 c(3) c(2) c(5) c(17) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(4) 0 0 0 0 0 0 0 0 0 0 a(6) a(16) 0 0 0 0 0 0 0 0 0 0 0 a(17) a(5) 0 a(14) a(8) a(4) a(12) 0 0 0 0 0 0 0 0 b(5) b(6) 0 0 b(12) 0 b(17) b(8) 0 0 b(4) 0 b(14) 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(12) 0 0 0 c(16) 0 0 c(14) 0 0 0 c(3) c(4) 0 0 0 c(5) 0 c(15) c(18) 0 c(17) c(8) 0 c(6) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(3) 0 0 0 0 0 0 0 0 0 0 a(7) a(12) 0 0 0 0 0 0 0 0 0 0 0 a(18) a(15) 0 a(5) a(17) a(3) a(6) 0 0 0 0 0 0 0 0 b(15) b(7) 0 0 b(6) 0 b(18) b(17) 0 0 b(3) b(4) b(5) b(8) b(14) 0 0 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(6) c(14) 0 0 c(12) 0 0 c(5) 0 0 0 0 c(3) 0 c(4) 0 c(15) 0 0 0 c(8) c(18) c(17) 0 c(7) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(2) 0 0 0 0 0 a(19) 0 0 a(20) 0 0 0 0 0 0 0 a(17) 0 0 0 0 0 0 0 0 0 b(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 0 b(19) 0 0 0 b(17) 0 0 0 c(17) 0 0 0 0 0 0 0 0 0 c(20) 0 0 c(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(2) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(4) a(6) 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 a(17) 0 a(8) 0 a(19) 0 0 0 0 0 a(16) 0 0 b(6) b(17) 0 0 0 0 0 0 0 0 0 b(19) 0 b(8) 0 0 0 b(12) b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 c(12) 0 c(8) 0 0 0 c(2) c(19) 0 0 0 c(17) 0 c(18) 0 0 0 0 c(6) 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(4) 0 0 0 a(5) a(17) 0 0 0 a(12) 0 0 0 0 0 a(16) 0 a(14) a(19) a(4) 0 0 0 0 a(8) 0 0 0 0 0 0 0 b(5) b(4) b(17) 0 0 b(8) b(12) b(19) 0 0 0 0 0 0 0 0 0 b(14) 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 c(16) 0 0 c(14) 0 0 0 0 0 c(13) 0 0 0 0 c(4) 0 c(3) c(2) 0 c(19) 0 c(5) c(17) 0 c(12) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(6) 0 a(12) 0 a(17) 0 0 0 0 0 0 0 0 0 b(6) 0 0 0 0 0 0 0 0 0 b(17) b(8) b(12) 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 0 0 c(12) 0 0 0 c(18) c(17) 0 c(8) 0 c(6) 0 c(7) 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(12) 0 0 0 0 0 0 0 0 0 0 0 0 a(8) a(14) 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 b(14) b(12) 0 0 b(16) 0 b(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 0 0 0 0 0 0 0 c(4) 0 0 0 c(3) c(14) 0 c(5) c(17) 0 c(8) 0 0 c(12) 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(7) 0 a(6) 0 a(18) 0 0 0 0 0 0 0 0 0 b(7) 0 0 0 0 0 0 0 b(8) 0 b(18) b(17) b(6) 0 b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(12) 0 0 0 0 0 c(6) 0 0 0 0 c(18) c(8) c(17) 0 c(7) 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(10) 0 0 0 c(14) 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(9) a(13) a(20) 0 0 0 0 0 0 0 0 a(11) 0 a(3) 0 0 a(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) b(13) b(20) b(3) b(11) b(2) 0 0 c(11) c(3) c(2) 0 0 0 0 c(13) 0 0 0 c(20) 0 c(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(13) a(20) 0 0 0 0 0 0 0 0 0 a(3) 0 a(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(13) b(20) 0 b(2) b(3) 0 0 0 c(3) c(2) 0 0 0 0 0 c(20) 0 0 0 0 0 c(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) 0 0 0 0 0 0 0 0 0 0 a(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 b(2) 0 0 0 c(2) 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 u(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) 0 0 0 0 a(14) a(8) 0 a(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(4) 0 0 b(8) 0 0 b(10) 0 b(14) 0 0 0 0 0 b(4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 c(16) c(4) 0 0 0 0 0 0 0 0 0 c(8) 0 0 0 0 0 0 0 0 c(1) c(9) 0 0 0 0 c(10) 0 0 c(14) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(4) 0 0 0 0 0 0 a(10) a(4) 0 0 0 a(14) a(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(19) 0 a(16) 0 b(10) 0 b(4) 0 b(8) 0 b(14) 0 0 0 b(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 c(16) 0 c(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(9) c(13) 0 0 0 c(10) c(4) 0 c(14) c(8) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(4) a(19) 0 0 a(14) a(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 b(4) 0 b(19) b(14) 0 0 b(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 0 0 c(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(13) c(20) 0 0 0 c(4) c(19) c(14) c(8) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(19) 0 0 0 a(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(16) 0 0 0 0 0 0 b(19) 0 0 b(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(16) 0 0 0 c(16) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 0 0 c(19) 0 c(8) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 a(9) 0 0 a(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) b(1) b(13) 0 b(20) c(1) c(9) c(13) c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(3) a(7) 0 0 0 0 0 0 0 0 0 0 0 0 a(6) 0 a(18) 0 a(17) 0 a(2) 0 0 0 0 0 a(12) 0 0 b(7) b(18) 0 0 0 0 0 0 0 0 0 b(2) b(19) b(17) 0 b(8) 0 b(6) b(12) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(8) 0 c(12) 0 c(6) 0 c(17) 0 0 0 0 c(2) 0 c(19) 0 c(18) 0 0 0 0 0 0 c(7) 0 0 0 0 0 0; 0 0 0 0 0 0 0 u(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(5) 0 0 0 a(7) 0 0 0 0 0 0 0 0 a(15) 0 a(11) 0 0 0 a(6) 0 0 0 0 b(7) 0 0 0 0 0 0 0 0 b(11) 0 0 0 b(1) 0 b(10) 0 b(15) 0 0 0 b(5) b(6) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(7) c(15) 0 0 0 c(1) 0 0 c(5) 0 0 c(10) 0 c(11) 0 0 0 0 0 0 c(6); 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) a(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(1) 0 0 a(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(1) 0 b(9) b(20) b(13) 0 c(1) c(9) c(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 0 0 0 0 0 0 0 a(1) a(9) a(13) a(20) 0 a(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) 0 0 a(4) b(20) 0 0 0 b(1) b(13) 0 b(9) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(10) 0 b(4) 0 b(19) 0 c(10) c(4) c(19) 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(1) c(9) c(13) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(9) a(13) a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(10) 0 a(4) 0 0 a(19) 0 0 0 0 b(9) b(20) 0 b(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(4) b(10) b(19) 0 0 c(10) c(4) c(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(9) c(13) c(20) 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(13) a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(4) 0 a(19) 0 0 0 0 0 0 0 b(13) 0 0 b(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(19) b(4) 0 0 0 c(4) c(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(13) c(20) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(19) 0 0 0 0 0 0 0 0 0 b(20) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(19) 0 0 0 c(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 u(2) 0 0 0 0 a(13) a(20) 0 0 0 0 0 0 a(11) a(3) a(2) 0 0 0 0 0 0 0 0 0 0 a(15) 0 a(9) a(18) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(9) b(13) b(20) 0 b(11) b(3) b(2) b(15) 0 b(18) 0 0 0 c(15) c(18) 0 0 0 0 c(3) c(9) c(13) 0 c(2) 0 c(11) 0 0 0 0 0 0 0 0 c(20) 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) 0 0 0 0 0 0 0 a(3) a(2) 0 0 0 0 0 0 0 0 0 a(15) 0 a(18) 0 a(13) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(13) b(20) 0 0 b(3) b(2) 0 b(18) b(15) 0 0 0 c(15) c(18) 0 0 0 0 0 c(2) c(13) c(20) 0 0 0 c(3) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(13) a(20) 0 0 a(4) a(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(14) 0 a(8) 0 0 0 0 b(13) 0 b(20) b(4) 0 0 b(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(8) b(14) 0 0 0 c(14) c(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(13) c(20) c(4) c(19) 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(20) 0 0 0 a(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 a(8) 0 0 0 0 0 0 b(20) 0 0 b(19) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 b(8) 0 0 0 c(8) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 c(20) 0 c(19) 0 0 0 0]';
end

function [H] = Hessian(f1coeff, f2coeff, f3coeff, s)
% the vector of monomials is
% m = [ const ; s1^2 * s2 ; s1 * s2 ; s1 * s3 ; s2 * s3 ; s2^2 * s3 ; s2^3 ; ... 
%       s1 * s3^2 ; s1 ; s3 ; s2 ; s2 * s3^2 ; s1^2 ; s3^2 ; s2^2 ; s3^3 ;   ...
%       s1 * s2 * s3 ; s1 * s2^2 ; s1^2 * s3 ; s1^3] 
%

% deriv of m w.r.t. s1
Hs1 = [0 ; 2 * s(1) * s(2) ; s(2) ; s(3) ; 0 ; 0 ; 0 ; ...
       s(3)^2 ; 1 ; 0 ; 0 ; 0 ; 2 * s(1) ; 0 ; 0 ; 0 ; ...
       s(2) * s(3) ; s(2)^2 ; 2*s(1)*s(3); 3 * s(1)^2];

% deriv of m w.r.t. s2
Hs2 = [0 ; s(1)^2 ; s(1) ; 0 ; s(3) ; 2 * s(2) * s(3) ; 3 * s(2)^2 ; ...
       0 ; 0 ; 0 ; 1 ; s(3)^2 ; 0 ; 0 ; 2 * s(2) ; 0 ; ...
       s(1) * s(3) ; s(1) * 2 * s(2) ; 0 ; 0];
   
% deriv of m w.r.t. s3
Hs3 = [0 ; 0 ; 0 ; s(1) ; s(2) ; s(2)^2 ; 0 ; ...
       s(1) * 2 * s(3) ; 0 ; 1 ; 0 ; s(2) * 2 * s(3) ; 0 ; 2 * s(3) ; 0 ; 3 * s(3)^2 ; ...
       s(1) * s(2) ; 0 ; s(1)^2 ; 0]; 

   
H = [ f1coeff' * Hs1 , f1coeff' * Hs2 , f1coeff' * Hs3;
      f2coeff' * Hs1 , f2coeff' * Hs2 , f2coeff' * Hs3;
      f3coeff' * Hs1 , f3coeff' * Hs2 , f3coeff' * Hs3];
end

%function [C] = cayley2rot(s)
%    C = ( (1-s'*s) * eye(3) + 2 * skewsymm(s) + 2 * s * s')' / ( 1 + s' * s);
%end